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Abstract 

Mass tensor molecular dynamics was first introduced by Bennett [J. Comput. Phys. 19, 267 
(1975)] for efficient sampling of phase space through the use of generalized atomic masses. Here, 
we show how to apply this method to ab initio molecular dynamics simulations with minimal 
computational overhead. Test calculations on liquid water show a threefold reduction in compu- 
tational effort without making the fixed geometry approximation. We also present a simple recipe 
for estimating the optimal atomic masses using only the first derivatives of the potential energy. 
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I. INTRODUCTION 



The efficiency of molecular dynamics (MD) simulations is governed by the choice of time 
step, which is usually taken 1/15 - 1/20 of the period of the fastest oscillation in the system. 
In liquid water, for instance, the O-H stretching mode has a period of 9-10 fs, thus limiting 
the time step to less than 1 fs. On the other hand, the length of a simulation must be 
much longer than structural relaxation times of the system, which range from picoseconds 
for liquid water to microseconds for biological systems [l|. Therefore, much effort has been 
made to reduce this gap in time scale f^yj], e.g., by eliminating the high frequency motion, 
or equivalently, by accelerating the low frequency motion. 

One of the most widely used technique in classical MD is constrained molecular dynamics 
(CMD) (4, 5], in which the bond lengths and angles may be fixed to their equilibrium values. 
CMD allows us to increase the time steps by a factor of two to three, due to the absence 
of fast intramolecular vibrations. CMD is equally valid in ab initio molecular dynamics 

n □ 

(AIMD) simulations [6|, |7[ which provide more accurate interatomic forces at the expense 



of much higher computational costs 



4l0j. We can also use CMD to 



arevent (undesirable) 



breaking of covalent bonds in the initial stages of the simulations 



llj . However, bond 



lexibility can play an important role if two or more phases coexist in the simulation cell 
I2I . I13I ]. Furthermore, overuse of bond angle constraints is prone to numerical artifacts j^J. 

A more recent development, which is free from these drawbacks, is the multiple time 
step (MTS) algorithms like Verlet-I [14] and r-RESPA [15]. In MTS, one can use large 
time steps for computationally expensive long-range interactions, while inexpensive short- 
range ones are integrated with small time steps. MTS is superior to CMD in that no fixed 
geometry approximation is required, while a similar or higher gain in performance can be 
achieved. Unfortunately, interatomic forces from ab initio calculations cannot be divided into 
short- and long-r ang e components exactly. Therefore, with the exception of approximate 
implementations flM, MTS is not generally used for AIMD. 

In the present paper, we investigate an alternative approach which takes advantage of the 
extra degrees of freedom associated with the choice of atomic masses. A simple approach 
along these lines is to rescale the atomic masses appropriately, which is particularly effective 
in systems containing hydrogen 18M2l|. as well as for adiabatic free energy calculations 
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23]. It is also a common practice to replace hydrogen by deuterium in AIMD study of 



liquid water |24|, thereby increasing the time step by a factor of \/2. 

Whereas the mass scaling method alleviates the time-scale problem to some extent, it does 
not fully eliminate the strong anisotropy of the phase space caused by covalent bonds. This 
problem can be overcome by using a more general formulation termed mass tensor molecular 
dynamics (MTMD), which was first introduced by Bennett more than three decades ago 25]. 
MTMD enables us to make the phase space nearly isotropic by using a nondiagonal position- 
dependent mass tensor, thus leading to enhanced sampling. Bennett used a constant, but 
nondiagonal mass tensor in his calculations on a simple polymer chain, and achieved a five- 
;o tenfold reduction in computational effort. The staging method used in path integral MD 
9] can also be viewed as a variant of MTMD. Melchionna has recently developed a similar 



method using internal coordinates in a series of papers [26l-l28l|. In the present work, we 
show how to apply MTMD to AIMD simulations with minimal computational overhead. 
We also demonstrate the effectiveness of our approach in Born-Oppenheimer MD simulation 
of liquid water. 



II. THEORY 



The classical Hamiltonian for a system of iV atoms is given by 



(1) 



where m, q and p are vectors of dm™ 3iV, representing atomic masses, positions and 
momenta. U(q) denotes either the classical potential energy function [29], or the Kohn-Sham 



total energy 
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31[. In what follows, we will focus on the latter prob 



The generalized Hamiltonian introduced by Bennett is given by [25 



em. 



# (q, P) = ^ X> Mr.\q) Pj + C/(q) 



(2) 



where M is a positive definite, symmetric matrix of dimension 3N which depends on the 
atomic positions q. Let us assume that p(q) is the canonical probability distribution in 
configuration space: 

/ dp exp(-pH) 



p(q) 



f dqdp exp(— j3H) 



(3) 



where = l/k-oT is the inverse temperature. Then, the ensemble average of any operator 
A(q) is given by 

{A)= f A(q)p(q)dq. (4) 

3211. we obtain 



If Eq.flSJ) is substituted into Eq.(j3J), and integration over p is performed 

= |M(q)|^exp(-/^(q)) 
Jdq|M(q)|lexp(-^(q))' 

Therefore, as long as det(M) = |M(q)| is independent of q, p(q) (and thus (A)) agrees with 
that of the original Hamiltonian. The same conclusion holds true for the microcanonical 
ensemble 



25) 



There are numerous ways to choose M that satisfies the above conditions, but a reasonable 
choice is 

M^-Hxc(q), (6) 
where the Hessian matrix H is given by 

«« = » (T) 
oqidqj 

and c(q) is a scalar function to compensate for the q-dependence of [H\. If Z7(q) is a convex 
quadratic function of q, Eq. ([6]) is optimal in the sense that all normal modes have the same 
frequency. In practice, it is prohibitively expensive to calculate the exact values of 7i%j at 
each step of AIMD. Instead, we introduce a simple harmonic approximation to U(q), given 
by 

angles bonds 

where only covalent bonds contribute to the sum. Then, we can define 

which is a sparse, symmetric, and indefinite matrix. The simplest way to construct a positive 
definite matrix similar to V is to define 

M = (V + el) x c(q) (10) 



using a small positive constant e. In our experience, however, the value of e tends to 
be relatively large in finite temperature simulations, which has a negative impact on the 
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performance of MTMD. Therefore, we adopt a more robust definition based on the eigenvalue 

n 

decomposition of V [25] . V is first diagonalized by an orthogonal matrix P as 



V = P T ■ A ■ P, (11) 

where A is a diagonal matrix given by Ay = Then, we introduce 

M = P T ■ /(A) • P, (12) 

where /(A) is a filtering function given by 

/(A) = (A 6 + 6 6 )i (13) 

At variance with Eq.f llOp . Mq is a positive definite matrix for any nonzero e. Finally, M is 
given by 

M = c M , (14) 

where 

c -(m)" (16) 

and a is an arbitrary normalization factor which does not depend on q. By definition, 
det(M) = a holds for any q, irrespective of the choice of Vr and e. 

The equations of motion derived from the generalized Hamiltonian of Eq.(J2D are given by 

^r&^'w d6) 

and 

di <9g fc 2 ^ % J <9g fc 



Numerica 



28|,l33r|35| 



integration of these equations is performed with the generalized leapfrog algorithm 



c +1 = # + 1 5>" +l {( M ^) n+1 + W 1 } > (19) 
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where h is the time step and the superscript denotes the time-step number. While Eq. (l20|) 
is an explicit formula, Eqs. ffTB"]) and ffl~9]) are implicit, i.e., these equations must be solved 
iteratively. Note, however, that we need to update only classical variables during the itera- 
tions, which is orders of magnitude faster than the evaluation of quantum mechanical forces 
(—dll/dq). Therefore, the computational cost of a single MTMD step is comparable to that 
of a conventional AIMD. Eqs. (jT8l - [20]) also preserve the symplectic property, if the implicit 



ones are solved with sufficient accuracy 



351 ] . When M is indepedent of q, the generalized 



leapfrog algorithm reduces to the velocity- Verlet method 



III. COMPUTATIONAL DETAILS 



We carried out two AIMD simulations of liquid water to evaluate the performance of 
MTMD in real applications. The reference simulation (hereafter denoted by REF) was per- 
formed with m H = 1.00794 and m = 15.9994, while the MTMD simulation was performed 
with q-dependent masses, as will be explained below. 

Atomic forces were calculated within the density functional theory 3711381. We used the 



generalized gradient approximation in the Perdew-Bur 



norm-conserving pseudopotentials were employed 40 



te-Ernzerhof form [39j . The separable 



411 ] . and only the T point was used 



to sample the Bril 



basis functions 



42 



ouin zone. The electronic orbitals were expanded by the finite-element 



431 ] with an average cutoff energy of 58 Ry, while t 



re resolution was 



approximately doubled near the oxygen atoms by adaptation of the grids 



44J. Liquid water 



at an elevated temperature was modeled by 64 molecules in a cubic supercell of side 13.92 
A 45], which was chosen to minimize the effect of nonergodic behavior observed at low tem- 



peratures 



46]. The equations of motion for the atoms were integrated with the generalized 



leapfrog algorithm using a time step of 20 a.u. (0.484 fs). Eqs.f lT8|) and f fT9|) were iterated 
10 times to achieve full convergence in the MTMD run. Initial atomic velocities were chosen 
so that the total energies in the two simulations coincide. After equilibration, production 
runs of 15 ps were carried out in the microcanonical ensemble. The electronic states were 
quenched to the Born-Oppenheimer surface at each MD step with the limited- memory vari- 



ant 
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48| of the quasi-Newton method (49 1. The initial guesses were extrapolated from 



previous MD steps [50 ] . 

Atomic masses used in the MTMD run were determined in the following manner. V% of 



Eq.flH]) was first decomposed into a sum of contributions from each water molecule, 



Vfe = £>/«• (21) 
The potential energy function t> M for molecule can be written as 

v, = y {K, - r ) 2 + (r 2 ,„ - r ) 2 } + - £ ) 2 , (22) 

where r x = r(OHi), r 2 = r(OH 2 ), and 6 = ZHiOH 2 . The parameters k±, A; 2 ,r , and Oq were 
determined by the force matching method |51] using the equilibration part of the trajectories, 
as listed in Table [B Under this definition, M is a block diagonal matrix, with each block 
being a square matrix of order 9. Therefore, inversion and diagonalization of M can be 
carried out at negligible cost. In terms of the efficiency of phase space sampling, the value 
of e introduced in Eq. ( 113ft should be as small as possible, while too small a value may lead 
to instability. After some trial and error, e = 0.1 Ha/Bohr 2 proved to be a good compromise 
in the present system. The normalization factor a was adjusted so that the period of the 
fastest oscillation in MTMD agrees with that in REF, to make a fair comparison between 
the two simulations. 



IV. RESULTS 

A. Numerical accuracy 

We first show the time evolution of total energy and potential energy during the MTMD 
run in Fig{T](a), which reflects the accuracy of numerical integration. Conservation of the 
total energy in REF is also satisfactory. Average temperatures are 415.74 K and 415.42 K 
for the REF and MTMD run, respectively. Probability distributions of potential energies 
are also compared in Fig{TJ(b). These results suggest that, within statistical errors, the two 
simulations are sampling the same region of phase space. We have also found that 20 - 25 % 
more computational effort is required for the calculation of interatomic forces in going from 
REF to MTMD. This is because the atoms move a longer distance at each MTMD step, as 
will be shown below, which makes the extrapolation of initial guesses less effective [50j. A 

n 

similar effect was observed in our previous work on CMD [6j. 
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B. Structural properties 



Average geometry of each molecule and the radial distribution functions (RDFs) have 
been calculated from the trajectories, as shown in Table [TT] and Figf3[ respectively. All 
results are in excellent agreement with each other, as expected from the theoretical analysis. 
We also compare the convergence of oxygen-oxygen RDFs (g 00 ) with respect to simulation 
length. The residual error in g OQ is defined by 

/*f max 

R (t) = \9oo(r,t) - g 00 (r,t max )\ 2 dr, (23) 

Jo 

where we set r max = 7 A and t max = 15 ps, and g o( r ,t) denotes the oxygen-oxygen RDF, 
extracted from the trajectory in the range of [0,t]. As shown in FigJS^a), the error for 
MTMD decays more than three times faster. The rapid convergence of MTMD is more 
clearly seen in Figj3]^b), where the RDFs at t = 1 ps are compared with that of t = t max . 
The residual errors for oxygen-hydrogen and hydrogen-hydrogen RDFs behave similarly, but 
are much smaller because there are twice as many hydrogen atoms in the simulation cell. 



C. Dynamical properties 

The atomic velocities from the simulations have been used to calculate the vibrational 
spectra shown in FigJH The first peak at 3400-3700 cm -1 in REF corresponds to the O-H 
stretching mode, while the second peak at 1620-1630 cm -1 is assigned to the H-O-H bending 
mode. The low frequency region below 1000 cm -1 corresponds to the intermolecular modes. 
When going from REF to MTMD, the two peaks corresponding to the intramolecular modes 
merge into a single broad peak between 3200 - 3700 cm -1 . Furthermore, the intermolecular 
modes now extend up to 3000 cm -1 . Therefore, the gaps in the original spectra, which 
reduce the efficiency of simulations, are absent in the MTMD results. Note also that the 
highest frequency observed in MTMD is in good agreement with that in REF, which confirms 
that our choice of a is appropriate. This also means that MTMD and REF have the same 
accuracy, in the sense that the time step ~ 1/20 of the period of the fastest oscillation 
(~ 9 fs) in both simulations. Finally, we compare the mean square displacements of oxygen 
atoms from the two simulations in FigEJ In agreement with the convergence rate of RDFs, 
the mean square displacement from MTMD is found to grow m 3.5 times faster. If the 
aforementioned increase in computational cost is taken into account, the use of MTMD 
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results in a net gain of 2.8 - 2.9 times. This gain is competitive with that of using a rigid 
water model |7], even though no fixed geometry approximation is required in MTMD. 



V. DISCUSSION 



While ensemble averages of structural properties remain unchanged in MTMD, the same 
does not hold for dynamical properties like vibrational spectra and self-diffusion coefficients. 
This is because the trajectories satisfying the generalized equations of motion (Eqs. (ll6|17p ) 
are unphysical, which can be viewed as a trade-off for enhanced sampling of the phase space. 
Therefore, MTMD should be used only in the equilibration phase, if dynamical properties 
are required. Nevertheless, there are several possible ways to recover the vibrational spectra 



of the original system within the harmonic approximation 
rely on the covariance matrix of atomic fluctuations, 



52M55| . Most of these methods 



a%j = <(?<- (Qi))(Qj -(Oj))), 



(24) 



which has the form of an ensemble average. Therefore, we can use the trajectory of MTMD 
or even Monte Carlo simulations to calculate Eq. (|24p . In disordered systems, however, care 
must be taken in the choice of reference frame to take into account the effect of translation 
and rotation of each molecule. It is also worth noting that a novel algorithm based on 
the stochastic path-integral formalism 56j has recently been proposed, which allows us to 
recover the time correlation functions of the original system from the trajectory on a mod- 
ified potential energy surface. A more empirical approach for estimating the self-diffusion 
coefficient from RDFs is also known for simple liquids 57j|. The validity of these algorithms 
will be discussed in more detail elsewhere. 

As already mentioned in the Introduction, MTMD generally performs better than mass 
scaling, because the knowledge of the molecules can be fully exploited in MTMD. Never- 
theless, mass scalin g is still an attractive option, if chemical reactions can occur during the 



simulations [ll|, |58|, |59|. Moreover, at variance with MTMD, no programming effort is re- 
quired. To the best of our knowledge, relatively little effort has been made 60] to optimize 
the atomic masses theoretically for finite temperature simulations, although several empir- 



ical studies exist 



19 



23]. To this end, we first introduce the expression for the thermal 



average of the Hessian |36|, |55j : 



CHij) 



d 2 U 
dq i dq j 



(9U\ fdU 



k B T \\dqij \dq 



(25) 



which was originally intended for vibrational analysis 55[. Here we propose to use Eq. (j25p 
to estimate the optimal masses m* heory by 



m 



theory 



oc (Ha) 




(26) 



where we take the average over all directions and atoms of the same element. Then, these 
masses may be used in conjunction with the original Hamiltonian of Eq.([T]). Eq. (!26|) may 
be viewed as a special case of Eq.(j6]). In the case of liquid water, we obtain a ratio of 2.0 
for mo heory /mH heory , which is in close agreement with the empirical results of Feenstra et 



al. 



19[. Theoretical prediction will play an important role in more complicated systems 



containing three or more elements where empirical optimization is prohibitive. Preliminary 
studies on the lithium borohydride (LiBH 4 ) system 61] using the optimal masses are showing 
promising results. We also note that atoms of the same element may have different masses. 
For instance, hydrogen atoms forming O-H bonds may be made heavier than those forming 
C-H bonds, if these bonds remain stable throughout the simulations. 



VI. CONCLUSION 



he importance of the choice of Hessian in ah initio geometry optimization is well known 



62 



65]. As we have shown in the present study, Hessian plays an equally important role 



in finite temperature simulations. MTMD will be the method of choice for nonreactive 



systems, if the primary interest of the simulations is the equilibrium properties 66] . On the 



other hand, the mass scaling method is applicable to a wider class of problems, and will 
be particularly effective for systems with large differences in atomic masses, e.g., water/Pt 
interface jfijj]. These approaches may also prove useful for global optimization problems 



681 ] . It would also be interesting to combine MTMD with rare event methods such as 
hyperdynamics [69'j and metadynamics (3]. Since these methods do not involve kinetic part 
of the Hamiltonian, the implementation would be straightforward. 
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FIG. 1. Numerical accuracy of test calculations on liquid water, (a) Time evolution of total energy 
and potential energy in MTMD run. (b) Probability distributions of potential energies from REF 
and MTMD runs. 
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FIG. 3. Convergence of oxygen-oxygen RDFs (goo( r ,t)) with respect to simulation length t. (a) 
Residual error R{t) as a function of time. The scaled MTMD lines denote R(t/3.5). (b) Comparison 
of g o(r, t) at t = 1 ps and t = t max . 
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FIG. 4. Vibrational spectra of liquid water from REF and MTMD. 
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FIG. 5. Mean square displacements of oxygen atoms from REF and MTMD. No average over time 
is taken. 
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TABLE I. Parameters of the force field for water molecules determined by the force matching 
method [511 ]. 



ki (kcal/mol-A 2 ) fc 2 (kcal/mol-rad 2 ) ro (A) 9q (deg) 



846.87 



74.33 



0.9846 105.01 



TABLE II. Average structure of water molecules in liquid phase from the simulations. 

r(OH) (A) ZHiOH 2 (deg) 
REF 0.9856 ± 0.0308 104.80 ± 5.72 
MTMD 0.9857 ± 0.0312 104.84 ± 5.75 
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